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The unstructured tetrahedral grid cell-centered finite volume flow solver USM3D has 
been recently extended to handle mixed element grids composed of hexahedral, prismatic, 
pyramidal, and tetrahedral cells. Presently, two turbulence models, namely, baseline 
Spalart-Allmaras (SA) and Menter Shear Stress Transport (SST), support mixed element 
grids. This paper provides an overview of the various numerical discretization options 
available in the newly enhanced USM3D. Using the SA model, the flow solver extensions are 
verified on three two-dimensional test cases available on the Turbulence Modeling Resource 
website at the NASA Langley Research Center. The test cases are zero pressure gradient flat 
plate, planar shear, and bump-in-channel. The effect of cell topologies on the flow solution is 
also investigated using the planar shear case. Finally, the assessment of various cell and face 
gradient options is performed on the zero pressure gradient flat plate case. 
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Nomenclature 

two-dimensional 
speed of sound 
boundary condition 
total drag coefficient 
pressure drag coefficient 
viscous drag coefficient 
local skin friction coefficient 
lift coefficient 
pressure coefficient 
grid spacing 
characteristic length 
Mach number 

number of control volumes or cells in a grid 
static pressure 
Reynolds number 
root mean square 

Spalart-Allmaras turbulence model 
Shear Stress Transport turbulence model 
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static temperature 
velocity magnitude 
streamwise velocity 

inner wall variable for streamwise velocity 

streamwise coordinate direction 

inner wall variable for distance 

coordinate direction normal to the aerodynamic body 

eddy viscosity 

laminar viscosity 


Subscripts 

ref = reference condition 

t = stagnation condition 

oo = freestream condition 


I. Introduction 

The NASA Tetrahedral Unstructured Software System (TetrUSS) [1, 2] was developed during the 1990s to 
provide a rapid aerodynamic analysis and design capability to applied aerodynamicists. The system is comprised of 
loosely integrated, user-friendly software that enables the application of advanced Euler and Navier- Stokes 
tetrahedral finite volume technology to complex aerodynamic problems. The system consists of component software 
for setting up geometric surface definitions (GridTool) [3], generating tetrahedral grids (VGRID) [4-6], computing 
Euler and Navier-Stokes flow solutions (USM3D) [7-12], and extracting meaningful information from analysis of 
results (SimpleView). The system is also coupled with the design tool CDISC [13]. Significant extensions in 
VGRID and USM3D capability are described in Ref. [2]. 

The USM3D flow solver has been widely used over the last two decades within NASA [14, 15], U.S. 
government [16] and industry [17] as a workhorse for aerodynamic analysis of complex configurations. A recent 
example is its role as a lead CFD flow solver within the Ares project of the NASA Constellation Program, 
furnishing over 7,000 solutions for the vehicles up through supersonic speeds [14,15]. USM3D was utilized in the 
aerodynamic development of the Crew Exploration Vehicle and its Launch Abort System configuration, providing 
over 1,000 solutions [2]. USM3D is also being used extensively in the NASA Aviation Safety Program to generate 
dynamic stability and control databases for civil transports for stall/post-stall flight regimes [18]. 

Some of the speed of USM3D computations is attributed to an analytical formulation of the spatial differencing 
stencil that exploits invariant features of tetrahedra [8,19], which renders it unnecessary to explicitly compute and 
store the cell gradients. This advantage is complimented by the VGRID tetrahedral grid generator, which generates 
organized thin layers of tetrahedral cells in the near-wall boundary layer region using Advancing Layer Method 
(ALM) [5], as illustrated in Fig. 1. Unfortunately, the diagonal faces and skewed topology of the highly stretched 
tetrahedra result in a severely distorted spatial differencing stencil for computing fluxes, which increases the 
numerical dissipation within the near-wall flow solution. Other cell topologies, such as hexahedra and prisms do not 
manifest this property. 

Since many recent applications have ventured into predominately separated flow regimes, attention is warranted 
on the sensitivity of smooth surface flow separation, such as a wing leading edge near stall, to near-wall cell 
topology. A similar case can be made for computing heat transfer on high-speed bodies. The accuracy of USM3D 
simulations for such flows can be improved by using more flow-aligned hexahedral or prismatic cells in the 
boundary layer that transition to tetrahedral cells away from the surfaces. 

Work is underway to expand the available cell topologies within USM3D to include hexahedra, prisms, and 
transitional pyramids. This paper describes those extensions and presents results from an initial code verification 
study. This study utilizes three two-dimensional test cases available on the Turbulence Modeling Resource website 
at the NASA Langley Research Center; the zero pressure gradient flat plate, planar shear, and bump-in-channel. 
Code verification is first performed on hexahedral grids for direct comparison to solutions from a well-established 
structured-grid code. The effect of cell topologies on the flow solution is assessed using a planar shear case. Finally, 
the characterization of the various cell and face gradient options is performed on the zero pressure gradient flat plate 
case. 
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II. Overview of USM3D 


The current USM3D [7, 8] code is a parallelized tetrahedral cell-centered, finite volume Navier-Stokes flow 
solver. The term “cell centered” means that the flow variables are solved at the centroid of each tetrahedral cell. 
In viscid flux quantities are computed across each tetrahedral cell face using various upwind schemes. Spatial 
discretization is accomplished by a novel reconstruction process, based on an analytical formulation for computing 
solution gradients within tetrahedral cells. The solution is advanced in time by a 2 nd -order Newton time step scheme 
[9], or to a steady-state condition by an implicit backward-Euler scheme. Several turbulence models are available: 
the Spalart-Allmaras (SA) one-equation model, the two-equation k-s turbulence model, the Menter Shear Stress 
Transport (SST) two-equation model, and the nonlinear Algebraic Reynolds Stress Models (ARSM) of Girimaji and 
Shih/Zhu/Lumley [8, 10, 11]. Detached Eddy Simulation (DES) has been implemented in all of the turbulence 
models. A capability to trip the flow at specified locations on aerodynamic surfaces has been implemented for the k- 
8 turbulence model, but fully turbulent flow is assumed for the results to follow. USM3D has capabilities for 
dynamic grid motion and overset grids [12]. 

III. Mixed Element Extensions to USM3D 

For second-order spatial accuracy in a cell-centered finite volume Navier-Stokes flow solver, gradients of 
solution variables are needed within a cell as well as on the faces of a cell. The gradients are used for the 
computation of convective and diffusive fluxes for the mean flow and turbulence model equations, as well as for the 
turbulence model source terms. Several changes have been incorporated in the tetrahedral grid version of USM3D to 
extend it for mixed element grids. Solution reconstruction based on an analytical formulation [19] that is applicable 
only for tetrahedral cells has been substituted by an explicit Taylor series expansion of primitive variables within 
each cell. Consequently, the cell gradient of a solution variable is explicitly evaluated and stored in solver memory. 
The extended flow solver allows flexibility to compute cell and face gradients using any of the several methods 
described below. In this section, key elements of solution methodology will be summarized. 

A. Geometric Quantities 

Cell center coordinates are calculated from arithmetic averaging of the coordinates of nodes that define a cell. 
Similarly, face center coordinates are calculated from arithmetic averaging of the coordinates of face nodes. 
Directed area of a face is calculated from the directed areas of its sub-triangles formed using the face edges and the 
face center. Cell volume is calculated by summing up the volumes of its constituent tetrahedral cells. A constituent 
tetrahedron within a cell is constructed from a sub-triangle of a cell face and the center of the cell. 

B. Cell Gradient 

Several options have been implemented for the computation of a solution variable gradient within a cell. The 
options include Green-Gauss integration, unweighted or weighted least squares procedure, and Mitchell stencil [8]. 

In the Green-Gauss integration, solution variables at the center of a cell face are calculated based on the 
averaging of the solution variables at the face nodes. Solution variables at a node are interpolated from the cells 
surrounding the node using a weighted averaging procedure [19]. The weights in the averaging procedures can 
optionally be limited between 0 and 2 to avoid the violation of the principle of positivity and maintain solution 
stability. Figure 2 shows a schematic representation of an unstructured grid. With reference to Fig. 2, the Green- 
Gauss stencil for cell A consists of all the cells shown in the figure (cells A, 1-13). This stencil is called a fully 
augmented stencil. 

For a least squares procedure, cell gradient of a solution variable can be evaluated using stencils of differing size. 
For example, a nearest-neighbor stencil can be used that includes only the face-neighbor cells (cells A, 1, 2, 3). A 
nearest-neighbor stencil for a boundary cell comprises of ghost cells associated with all of its boundary faces. It is a 
very compact and computationally efficient stencil. However, for tetrahedral grids, a solution based on this stencil 
quickly becomes unstable. On the other hand, a fully augmented stencil consisting of all the cells surrounding each 
node of a cell can also be used in the least squares procedure. This stencil does not include ghost cells for a 
boundary cell. The stencil is stable for any cell topology, however it is computationally inefficient. An optimized 
smart augmentation [20] stencil that is stable yet efficient is also implemented, which is composed of cells A, 1, 2, 3, 
5, 8, and 12, in our example. An option to use fully augmented stencil for the boundary cells, and either nearest 
neighbor, fully augmented, or smartly augmented stencil for the interior cells is available as well. For the weighted 
least squares, contributions to a solution variable gradient from any cell within the stencil of a specific cell is 
weighted by the inverse of the distance between the two cell centers under consideration. 
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C. Face Gradient 

Face gradient of a solution variable can be evaluated from either Mitchell’s stencil [8] or an average of the 
gradients within two cells sharing the face. For the latter case, the cell-averaged gradient is suitably augmented using 
a face-normal augmentation [21], an edge-normal method [22], or a face-tangent [23] method. Various methods are 
illustrated next using a sketch below for a 2D grid. 

The sketch shows a grid face defined by two nodes pi and p2. The face is shared by cells C and D. A position 
vector f points from the center of cell C to the center of cell D. The unit normal to the face is represented by h. The 
face gradient computation procedure based on the Mitchell’s stencil is shown in Eq. (1) below. The Eqs. (2-4) 
illustrate face gradient based on the face-normal, edge-normal, and face-tangent methods, respectively. 
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r is the unit vector along the position vector r , such that r = r/|r|. In the Eqs. (1-4), V0y =<P x i +<P r J is the 
gradient of a scalar function <j> . ^(j> aV g represents the average of the gradients in cell C and D; 
'y&avg = (^0C + ^0£>)/2- 


D. Turbulence Models 

The Spalart-Allmaras [24] and the Menter Shear Stress Transport [25] turbulence models have been 
implemented to support mixed element grids. Various researchers have attempted several revisions to these two 
models. The standard forms of these models as well as some of the noteworthy revisions to them have been 
comprehensively documented in the NASA Langley Research Center Turbulence Model Resource website [26]. 
Following the nomenclature used in the website, “standard” Spalart-Allmaras One-Equation Model (SA), 
“Standard” Menter SST Two-Equation Model (SST), Menter SST Two-Equation Model with Vorticity Source Term 
(SST-V), and Menter 2003 SST Two-Equation Model (SST-2003) are presently available in the extended version of 
USM3D. 


IV. Verification Results 

USM3D extensions for mixed element grids are verified using three verification cases available on the 
Turbulence Model Resource website [26]. All the cases are quasi two-dimensional (2D), having one hexahedral cell 
in the spanwise direction. The cases are zero pressure gradient flat plate, bump-in-channel, and planar shear. For 
each case, a series of five nested structured grids in PLOT3D and CGNS formats are available on the Turbulence 
Model Resource website [26]. A utility is used to convert structured PLOT3D grids into unstructured hexahedral 
grids and write them in the VGRID .gridu format. USM3D solutions are computed based on the unstructured 
hexahedral grids. Grid convergence of certain local as well as integrated quantities is monitored. Additional USM3D 
verification is achieved by code-to-code comparison. For this purpose, the finest grid USM3D and CFL3D [27] 
solutions are compared. The CFL3D results shown in this paper are obtained from the Turbulence Model Resource 
website. Verification results are presented in the sub- sections A, B, and C below. 

Grid topology effects on a solution are studied using a 2D planar shear case in sub-section D. For this purpose, 
grid convergence of various quantities from the hexahedral, prismatic, and tetrahedral grid solutions is presented. 
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Each of these grid types is generated from the same underlying structured grids. In sub-section E, three gradient 
computation options are assessed using the 2D zero pressure gradient flat plate case. 

USM3D computations reported in the present work treat turbulence model convective terms with first-order 
spatial accuracy. In viscid flux quantities are computed using Roe’s flux-difference- splitting (FDS) scheme. Steady- 
state solutions are obtained using local time stepping and implicit defect-correction scheme. USM3D solutions for 
all three verification cases have been computed using the “standard” Spalart-Allmaras One-Equation Model (SA) as 
well as the Menter SST Two-Equation Model with Vorticity Source Term (SST-V). However, for the sake of 
brevity, only the results from the SA model will be presented in the present paper. 

USM3D results in sub-sections A-D are based on the “GGM” option that uses Green-Gauss integration 
procedure for the computation of gradients within a cell and the Mitchell’s method for the computation of gradients 
at a face. In the GGM method, cell gradients are used for the evaluation of mean flow convective terms whereas the 
face gradients are used for the mean flow and turbulence model diffusive terms as well as the turbulence model 
source terms. 

A. 2D Zero Pressure Gradient Flat Plate 

For this case, five progressively finer unstructured hexahedral grids have been generated that are based on the 
underlying structured PLOT3D grids. The specific nested structured grids [26], denoted by the number of grid points 
in the span wise, streamwise, and body-normal directions, are 2x35x25, 2x69x49, 2x137x97, 2x273x193, and 
2x545x385. Solutions have been computed for the freestream Mach number of 0.2 and Reynolds number of 5xl0 6 
per unit length of the plate. The angle-of-attack is zero degree. Figure 3 presents a schematic view of the 
computational domain and the boundary conditions in the X-Z plane. Symmetry boundary condition is applied on 
the two spanwise planes. 

On each grid, a USM3D solution is computed starting from the freestream-based initial conditions. The solutions 
are computed using a progressively larger number of iterations to ensure accurate results on all grids. The number of 
iterations ranged from 25,000 for the coarsest grid to 150,000 for the finest grid. The convergence of the USM3D 
mean flow and SA model residual errors is presented in Fig. 4 for all five grids. It is evident from the figure that on 
all grids the rms of the mean flow residual errors is below 1.0e-15 and the SA model residual errors is below l.Oe- 
12 . 

USM3D verification results based on the SA model are presented next. Figure 5 shows grid convergence of the 
USM3D and CFF3D computed skin friction at a specific location on the plate as well as the drag coefficient 
integrated over the plate. It can be seen that as the grid is refined, USM3D and CFF3D results show closer 
agreement. Focal skin friction and the integrated drag from USM3D and CFF3D are almost identical on the finest 
grid. It should be noted that the coarsest grid USM3D and CFF3D results appear very different due to a narrow y- 
axis range. For example, on the coarsest grid USM3D and CFF3D drag values differ only by less than 1.2 count. 

In the next two figures the finest grid USM3D and CFF3D results are compared. The local skin friction 
variations along the plate are compared in Fig. 6, whereas, the profiles of streamwise velocity and eddy viscosity in 
the body-normal direction are presented in Fig. 7. The streamwise velocity profiles are presented in terms of the 
inner wall variables (u + and y + ). The profiles are presented at two streamwise stations on the plate (Fig. 7a). The 
eddy viscosity variations across the boundary layer are compared at one streamwise station on the plate (Fig. 7b). 
Figures 6 and 7 demonstrate almost identical results from USM3D and CFF3D. 

B. 2D Bump-in-Channel 

For this case, five progressively finer unstructured hexahedral grids have been generated that are based on the 
underlying structured PFOT3D grids. The specific nested structured grids [26], denoted by the number of grid points 
in the spanwise, streamwise, and body-normal directions, are 2x89x41, 2x177x81, 2x353x161, 2x705x321, and 
2x1409x641. Solutions have been computed for the freestream Mach number of 0.2, Reynolds number of 3xl0 6 per 
unit grid length, and angle-of-attack of 0 degree. Figure 8 presents a schematic view of the computational domain 
and the boundary conditions in the X-Z plane. Symmetry boundary condition is applied on the two spanwise planes. 

On each grid, USM3D solution is computed starting from the freestream-based initial conditions. Moving from 
the coarsest to the finest grid, the solutions required a progressively larger number of iterations to converge. USM3D 
is run for 40,000 iterations on the coarsest grid and for as many as 500,000 iterations on the finest grid. The 
convergence of the USM3D mean flow and SA model residual errors is presented in Fig. 9 for all five grids. The 
rms of the residual errors is reduced below 1.0e-12 for the mean flow equations and below 1.0e-6 for the SA model 
equation. On the finest two grids, the USM3D mean flow and SA model convergence is very slow, possibly due to 
the linear solver that is first-order accurate and uses point-implicit algorithm. This issue will be investigated in the 
follow-on work by using a consistent linearization and line-implicit scheme for the linear solver. 
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USM3D SA model verification results for this case are presented in Figs. 10-14. Grid convergence of the 
USM3D and CFL3D computed local skin friction coefficient at the bump apex and the integrated viscous drag 
coefficient is displayed in Fig. 10. Figure 1 1 presents the grid convergence of the total drag and lift coefficients from 
the USM3D and CFL3D solutions. Except for the viscous drag coefficient, all the quantities in Figs. 10 and 11 
exhibit similar convergence trends. The USM3D and CFL3D results on the coarsest grid differ the most. For 
example, the total drag from USM3D and CFL3D differs by about seven counts on the coarsest grid. However, as in 
the flat plate case, the skin friction and the integrated force coefficients converge to nearly identical values on the 
finest three levels of grid. 

Next, the finest grid USM3D and CFL3D solutions are compared. Figure 12 presents the streamwise variations 
of the skin friction and the pressure coefficient along the bump surface. In Fig. 13, streamwise velocity and eddy 
viscosity variations in the direction normal to the bump surface are shown at certain stations on the bump. The eddy 
viscosity field variations from the USM3D and CFL3D solutions are demonstrated in Fig. 14. For the sake of ready 
comparison, the contour plots in Fig. 14 use the identical range for both the solutions. The finest grid USM3D and 
CFL3D results as presented in Figs. 12-14 demonstrate an excellent agreement between the two flow solvers. 

C. 2D Planar Shear 

This case involves the growth of a free shear layer that develops when two different streams passing over a thin 
plate begin mixing. For this case, five progressively finer single-block unstructured hexahedral grids have been 
generated that are based on the underlying multi-block structured PLOT3D grids. Each structured grid consists of 
three blocks. In the description below, dimensions of all blocks will be listed parenthetically to represent a block- 
structured grid. The specific nested structured grids [26], denoted by the number of grid points in the span wise, 
streamwise, and body-normal directions, are (2x9x17, 2x9x17, 2x33x33), (2x17x33, 2x17x33, 2x65x65), (2x33x65, 
2x33x65, 2x129x129), (2x65x129, 2x65x129, 2x257x257), and (2x129x257, 2x129x257, 2x513x513). Figure 15 
presents a schematic view of the computational domain and the boundary conditions in the X-Z plane. Symmetry 
boundary condition is imposed on the two span wise planes. The Mach number of the smaller stream below the plate 
is 0.5, whereas, the Mach number of the larger stream above the plate is 0.25. Solutions have been computed for the 
Reynolds number of 50,000 based on unit grid length. 

For each grid, USM3D solution is computed starting from the freestream-based initial conditions. USM3D is run 
for 20,000 iterations each on the coarsest two levels of grids. For the next two finer levels of grid USM3D is run for 
30,000 and 50,000 iterations, respectively. On the finest grid, USM3D is run for 125,000 iterations. The 
convergence histories of the mean flow and SA model equation residual errors are similar to those for the flat plate 
case and therefore not shown here. It is noted that the rms of the USM3D residual errors on all grids is reduced 
below 1.0e-14 for the mean flow equations and below 1.0e-10 for the SA model equation. 

USM3D SA model verification results for this case are presented in Figs. 16-20. Grid convergence of the 
USM3D and CFF3D solutions is exhibited in Fig. 16 using two quantities, namely, the drag coefficient on the 
streams-dividing plate and the streamwise velocity at a specific location in the computational domain. The finest 
grid drag coefficient from USM3D and CFF3D differs by a little over Vi count. The local streamwise velocity from 
USM3D and CFF3D are seen to be converging to similar values. 

Similar to the previous two test cases, USM3D solution on the finest grid is assessed next using the 
corresponding CFF3D solution as the benchmark. Accordingly, the streamwise variations of streamwise velocity 
along the Z=0 line are compared in Fig. 17. The normal-direction variations of streamwise velocity and eddy 
viscosity at a specific streamwise station are presented in Fig. 18. Next, a qualitative comparison of the USM3D and 
CFF3D computed flowfields is made. For this purpose, eddy viscosity variations in a spanwise plane are shown in 
Figs. 19-20. The variations in the proximity of the streams-dividing plate are presented in Fig. 19, whereas the 
farfield variations are demonstrated in Fig. 20. Based on the Figs. 17-20, it is observed that the USM3D and CFF3D 
solutions on the finest grid are in excellent agreement. 

D. Effect of Grid Topology 

The grid topology effect on the computed solutions is assessed using the 2D planar shear case. Based on the 
family of five unstructured hexahedral grids used previously, three additional sets of grids are generated by 
subdividing hexahedral cells in various manners. Each of these sets consists of only one type of cells. Two sets 
consist of prismatic cells whereas the third one consists of tetrahedral cells. Two prismatic cell sets differ by the 
direction in which prism cells are aligned. In one set, prismatic cells are aligned in the direction normal to the 
streams-dividing plate. This set will be identified as “Normal prisms” in the discussion to follow. Another set, 
designated as “Sideways prisms”, comprises of prismatic cells that are aligned in the spanwise direction. The 
original hexahedral grid set will be called “Hexes”, whereas the one made up of tetrahedral cells will be named as 
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“Tets”. Using the finest grid in the respective sets, the characteristics of various sets are summarized in the Table 1 
below. 


Table 1. Statistics for the finest grid in various grid sets constructed for the planar shear case 


Grid Set Identifier 

Number of Nodes 

Number of Cells 

Number of 
Boundary Faces 

Number of Faces 

Hexes 

657,922 

327,680 

657,920 

1,312,000 

Normal Prisms 

657,922 

655,360 

659,456 

1,968,128 

Sideways Prisms 

657,922 

655,360 

1,313,280 

2,295,040 

Tets 

657,922 

1,966,080 

1,315,840 

4,590,080 


As in the previous cases, USM3D solution for each grid is computed starting from the freestream-based initial 
conditions. Following the grid notation system in Ref. [26], the five grids within each grid topology are identified by 
an integer number varying from 0 to 4, where 0 is used to denote the finest grid whereas 4 is used to denote the 
coarsest grid. Using this identification system the convergence characteristics of three grid levels within various grid 
topologies are outlined in Table 2. 


Table 2. Convergence characteristics of various grid sets constructed for the planar shear case 


Grid Set Identifier 

Grid Level 

Number iterations 

rms of mean flow 
residuals 

rms of turbulence 
model residuals 

Hexes 

4 

20,000 

1.84e-15 

5.96e-14 


2 

30,000 

5. Ole-16 

1.99e-14 


0 

125,000 

5.13e-16 

6.01e-10 

Normal Prisms 

4 

20,000 

4.20e-15 

7.87e-14 


2 

30,000 

5.12e-16 

1.12e-14 


0 

170,000 

1.27e-16 

1.24e-14 

Sideways Prisms 

4 

20,000 

2.03e-15 

3.21e-3 


2 

30,000 

5.00e-14 

6.46e-4 


0 

240,000 

2.67e-12 

6.87e-5 

Tets 

4 

20,000 

2.85e-15 

3.91e-4 


2 

60,000 

1.62e-14 

1.31e-4 


0 

310,000 

1.86e-10 

2.15e-5 


It is evident from Table 2 that the convergence behavior of “Normal Prisms” is similar to “Hexes” and “Tets” is 
similar to “Sideways Prisms”. The SA model convergence for the “Sideways Prisms” and “Tets” is unsatisfactory 
across all grid levels. Examination of the evolution of the SA model primary variable for the latter two grid 
topologies reveals that the variable rapidly approaches its artificially imposed minimum value of 1.0e-10 inside 
USM3D. The clipping of the SA model variable may be responsible for the stalled convergence. A separate study 
focusing on the implementation and evaluation of the “negative” SA model [28] is underway to examine a possible 
link between the artificial clipping of the SA model variable and a degraded convergence of the model. 

Grid convergence of USM3D solutions on various grid topologies is monitored in Fig. 21. For this purpose, the 
variations of drag coefficient on the streams-dividing plate and streamwise velocity at a specific location in the 
computational domain are plotted. It can be seen from Fig. 21 that with grid refinement all grid topologies converge 
to similar values. 

USM3D finest grid solutions from various grid topologies are examined next using the streamwise variations of 
streamwise velocity along the Z=0 line (Fig. 22) and the normal- direction variations of streamwise velocity and 
eddy viscosity at a specific streamwise station in the computational domain (Fig. 23). It is clear from Figs. 22-23 
that the finest grids solutions from various grid topologies are in close agreement. 

E. Assessment of Various Gradient Computation Options 

In this sub- section the 2D zero pressure gradient flat plate case is revisited to compare hexahedral grid solutions 
obtained by exercising various USM3D options for the computation of cell and face gradients. In addition to the 
USM3D “GGM” solutions reported in sub-section A, two additional USM3D solution sets have been computed for 
the same case. 
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Both of the additional sets use least squares procedure for the computation of gradients within a cell. One set, 
designated as “ULSQ”, uses unweighted least squares procedure whereas the “WLSQ” solution set uses the 
weighted least squares procedure. In both the sets, cell gradients are used for the evaluation of mean flow convective 
terms as well as the turbulence model source terms. In both sets, face gradients are computed from the average of 
respective cell gradients that are augmented by the “face-tangent” method. The mean flow and turbulence model 
diffusive terms are computed using these face gradients. 

On all five hexahedral grids, USM3D “ULSQ” and “WSLQ” solutions are computed using the same number of 
iterations as that used for the respective USM3D “GGM” solutions. The convergence of the mean flow and SA 
model in the solutions computed from the three options is very similar and therefore not shown here. For example, 
in the finest hexahedral grid USM3D “ULSQ”, “WLSQ”, and “GGM” solutions, the rms of the mean flow residuals 
ranged from 1.53e-15 to 1.54e-15 and the rms of the SA model residuals varied from 2.61e-13 to 2.77e-13. 

Figure 24 shows the grid convergence of skin friction at a specific location on the plate and the drag coefficient 
integrated over the plate, both of which are computed using various USM3D options described earlier. Local skin 
friction and integrated drag show close agreement among various USM3D solutions on the finest two grids. It can be 
seen from Fig. 24 that the local skin friction and the integrated drag from USM3D “WLSQ” and “GGM” solutions 
are in close agreement on all grids. However, the local skin friction and the integrated drag in the USM3D “ULSQ” 
solution are significantly lower on the coarsest grid. For example, on the coarsest grid, “ULSQ” drag is 
approximately three counts lower than the corresponding “WLSQ” and “GGM” values. 

In the next two figures the finest grid results from the three USM3D solutions are compared. Variations of local 
skin friction along the plate are plotted in Fig. 25, whereas, the profiles of streamwise velocity and eddy viscosity in 
the body-normal direction are presented in Fig. 26. The streamwise velocity profiles (Fig. 26a) at two streamwise 
stations are shown in terms of the inner wall variables (u + and y + ). The eddy viscosity variations (Fig. 26b) are 
compared across the boundary layer at one streamwise station on the plate. The results in Figs. 25 and 26 display 
almost identical agreement between the three USM3D solutions. 

V. Concluding Remarks 

The NASA USM3D tetrahedral grid cell-centered finite volume flow solver has been recently extended to 
accommodate mixed element grids composed of hexahedral, prismatic, pyramidal, and tetrahedral cells. The 
important work of code verification is well underway. The initial verification study that is the topic of this paper 
used three two-dimensional test cases available on the Turbulence Modeling Resource website at the NASA Langley 
Research Center; the zero pressure gradient flat plate, planar shear, and bump-in-channel. The approach was to first 
verify the implementation on hexahedral grids through direct comparison to solutions from the well-established 
structured-grid code CFL3D. Next the effect of cell topologies on the flow solution was assessed using the planar 
shear case. And finally, the various cell and face gradient reconstruction schemes were evaluated on the zero 
pressure gradient flat plate case using hexahedral grids. 

The integrity of coding implementations is verified through assessments of convergence for the flow and 
turbulence model equations, grid-convergence of both integrated and distributed pressure and skin-friction 
coefficients, and off-body comparisons of velocity and eddy viscosity distributions. The overall integrity of USM3D 
code modifications, which include a more general data structure, implementation of least-squares reconstruction, 
and the generalization of Green-Gauss scheme, is confirmed by near identical agreement with CFL3D solutions on 
hexahedral grids for all criteria. The only anomaly to surface was a stalled convergence on tetrahedral and 
“sideways” prisms for planar shear flow that may be caused by an artificial clipping inside the USM3D 
implementation of Spalart-Allmaras (SA) turbulence model. A separate study of a “negative” SA model is underway 
to investigate this problem. 

While the current verification study only included the SA model results, equally good correlations have been 
obtained with the Shear-Stress Transport (SST) model, but were not included for sake of brevity. These results will 
be subject of a future publication. 
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Figure 1. Thin-layered tetrahedral grid in the boundary layer region. 



Figure 2. A representative view of a 2D unstructured triangular grid. 
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P t /P ref = 1.02828 



Figure 3. Schematic view of computational domain and various boundary conditions for the 2D zero 
pressure gradient flat plate. 




Figure 4. USM3D convergence for the 2D zero pressure gradient flat plate on a sequence of five hexahedral 
grids. Mo 3=0.2, Re L =5xl0 6 (L=l). 
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Figure 5. Grid convergence of skin friction and drag coefficients for the 2D zero pressure gradient flat 
plate computed using SA model on a sequence of five hexahedral grids. M^O.2, Re L =5xl0 6 (L=l). 



Figure 6. Comparison of the streamwise variations of local skin friction computed using SA model on the 
finest hexahedral grid for the 2D zero pressure gradient flat plate. Moo=0.2, Re L =5xl0 6 (L=l). 
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Figure 7. Comparison of the normal direction variations of streamwise velocity and eddy viscosity computed 
using SA model on the finest hexahedral grid for the 2D zero pressure gradient flat plate. M^O.2, Re L =5xl0 6 
(L=l). 



x x 

Figure 8. Schematic view of computational domain and various boundary conditions for the 2D bump-in- 
channel. 
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Figure 9. USM3D convergence for the 2D bump-in-channel on a sequence of five hexahedral grids. M^O.2, 
Re L =3xl0 6 (L=l). 
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(b) integrated viscous drag coefficient 



Figure 10. Grid convergence of 2D bump-in-channel skin friction and viscous drag coefficients computed using 
SA model on a sequence of five hexahedral grids. M^O.2, Re L =3xl0 6 (L=l). 
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Figure 11. Grid convergence of 2D bump-in-channel total drag and lift coefficients computed using SA model on 
a sequence of five hexahedral grids. M^O.2, Re L =3xl0 6 (L=l). 




(a) skin friction coefficient (b) coefficient of surface pressure 

Figure 12. Comparison of the streamwise variations of skin friction and surface pressure coefficients computed 
using SA model on the finest hexahedral grid for the 2D bump-in-channel. Moo=0.2, Re L =3xl0 6 (L=l). 
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(a) streamwise velocity at two stream wise stations 


(b) eddy viscosity at X=0.97 


Figure 13. Comparison of the normal direction variations of the normalized streamwise velocity and eddy viscosity 
computed using SA model on the finest hexahedral grid for the 2D bump-in-channel. M^O.2, Re L =3xl0 6 (L=l). 



(a) USM3D 


(b) CFL3D 


Figure 14. Nearfield normalized eddy viscosity contours from USM3D and CFL3D computed using SA model on 
the finest hexahedral grid for the 2D bump-in-channel. Moo=0.2, Re L =3xl0 6 (L=l). 
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C D on dividing plate 



Figure 15. Schematic view of computational domain and various boundary conditions for the 2D planar 
shear case. 




Figure 16. Grid convergence of drag coefficient and local normalized streamwise velocity for the 2D planar shear 
case computed using SA model on a sequence of five hexahedral grids. Moo=0.5, Re L =50,000 (L=l). 
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Figure 17. Comparison of the normalized streamwise velocity along Z=0 station computed using SA model on 
the finest hexahedral grid for the 2D planar shear case. M oo =0.5, Re L =50,000 (L=l). 




(b) eddy viscosity at X=29.2468 


Figure 18. Comparison of the normal direction variations of the normalized streamwise velocity and eddy viscosity 
computed using SA model on the finest hexahedral grid for the 2D planar shear case. Moo=0.5, Re L =50,000 (L=l). 
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Figure 19. Nearfield normalized eddy viscosity contours from USM3D and CFL3D computed using SA model on 
the finest hexahedral grid for the 2D planar shear case. Moo=0.5, Re L =50,000 (L=l). 
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(a) USM3D (b) CFL3D 


Figure 20. Farfield normalized eddy viscosity contours from USM3D and CFL3D computed using SA model on the 
finest hexahedral grid for the 2D planar shear case. Moo=0.5, Re L =50,000 (L=l). 
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Figure 21. Grid convergence of drag coefficient and local normalized streamwise velocity for the 2D planar shear 
case computed using SA model on a sequence of five grids of four different topologies. Moo=0.5, Re L =50,000 (L=l). 



Figure 22. Comparison of the normalized streamwise velocity along Z=0 station computed using SA model on the 
finest grid from four different topologies for the 2D planar shear case. Moo=0.5, Re L =50,000 (L=l). 
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(a) streamwise velocity at X=29.2468 



Figure 23. Comparison of the normal direction variations of the normalized streamwise velocity and eddy viscosity 
computed using SA model on the finest grid from four different topologies for the 2D planar shear case. M^O.5, 
Re L =50,000 (L=l). 



(a) local skin friction at X=0.97 (b) integrated drag coefficient 

Figure 24. Grid convergence of skin friction and drag coefficients for the 2D zero pressure gradient flat plate 
computed using USM3D SA model and various gradient calculation procedures on a sequence of five hexahedral 
grids. Mco=0.2, Re L =5xl0 6 (L=l). 
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Figure 25. Comparison of the streamwise variations of local skin friction computed using USM3D SA model and 
various gradient calculation procedures on the finest hexahedral grid for the 2D zero pressure gradient flat plate. 
Mco=0.2, Re L =5xl0 6 (L=l). 




Figure 26. Comparison of the normal direction variations of streamwise velocity and eddy viscosity computed using 
USM3D SA model and various gradient calculation procedures on the finest hexahedral grid for the 2D zero pressure 
gradient flat plate. M^O.2, Re L =5xl0 6 (L=l). 
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